/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name :f_log.c
 *  Function  :
 *        f_log returns Logarithm value.
 *        This file is copied from math.lib of 
 *                  CC33 tool(CC33v40.exe + CCPA45.exe, 
 *                  math of ansilib33v25 )
 *
 *  original file's Revision  :
 *      2000/02/10    first release                 M.Igarashi
 *  Revision  :
 *      2003/04/10    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/
#include <f_math.h>
#include <f_smcvals.h>

#define	K0	-0.18327870372215593212e2f
#define	K1	 0.934639006428585382474e2f
#define	K2	-0.90174691662040536328986e2f
#define	S0	-0.20733487895513939345e2f
#define	S1	 0.61761065598471302843e2f
#define	S2	-0.45087345831020305748486e2f



// ALGORITHM                
//1. Argument Reduction: find k and xfrac such that 
//			x = 2^k * (xfrac), 
//	   where  sqrt(2)/2 < xfrac < sqrt(2) .
//
//2. Approximation of log(xfrac).
//	Let x = xfrac/(1+xfrac) ; based on log(xfrac) = log(1+x) - log(1-x)
//		 = 2s + 2/3 x**3 + 2/5 x**5 + .....,
//	     	 = 2s + x*R
//	2x = xfrac-1 - x*(xfrac-1) = (xfrac-1) - hfsq + s*hfsq, where hfsq = (xfrac-1)*(xfrac-1)/2.
//	In order to guarantee error in log below 1ulp, we compute log by
//		log(xfrac) = xfrac-1 - (hfsq - s*(hfsq+R)).
//	
//3. Finally,  log(x) = k*ln2 + log(xfrac).  
//
/* -- Format of floating point (single precision) ---
 * 
 * 31|30         23|22                0|
 *  -----------------------------------
 * | |   exponent  |   fraction        |
 *  -----------------------------------
 * |    8 bits          23 bits
 */



float f_log(float sfX){

	FLT_LNG       ckarg;       /* exchange the type to check the argument   */
	long lX;
	float	sfHfsq;
	float	sfXfrac,sfRet,sfTemp;
	int	iK;
	
	FLT_LNG		sfTemp2;
	
	F_GETW(lX,sfX);			// get x
	
/* check the argument */
	 ckarg._F = sfX;
	 F_CHECK_ARG( &ckarg );                   /* check NaN */
	                                          /* error: error code 33 domain error */
	if ( ckarg.st._LL == f_P_INF.st._LL ) {   /* x = Inf */
	    errno = ERANGE;                       /* 34: range error */
	    return( f_P_INF._F );
	} else if ( ( ckarg.st._LL == f_N_INF.st._LL ) || ( (lX&0x80000000) !=0 ) ) {  /* x < 0 */
	    errno = EDOM;                         /* 33: domain error */
	    return( f_NAN._F );
	} else if ( lX == 0 ) {
	    errno = ERANGE;                       /* 34: range error */
	    return( f_N_INF._F );
	} else if ( lX == 0x3f800000 ) {
		// for x = 1.0f
		// if below routine doesn't exist, log( 1 ) returns 0x34000000, 1.1920928955e-7.
		sfRet = 0.0f;
		return sfRet;
	}

	sfX = f_frexp( sfX, &iK );
	F_GETW(lX,sfX);			// get x
	sfTemp2.st._LL = 0x00000000;   //  set 0(float)
	if( (lX&0x80000000) !=0 || lX < 0x3f3504f3 ){    // means sfX < 1 / 2^(-0.5)

		// from comments, x should be below when calcuralting log function :
		//   1/sqrt( 2 ) < x < sqrt( 2 )
		//  but f_frexp function returns 0.5 - 0.9999999( not 1.0 ) .
		sfTemp2.st._LL = 0x3fb504f1;  //  nearly  sqrt( 2 )=0x3fb504f3, but smaller a little
		sfX *= sfTemp2._F;
		sfTemp2._F = 0.34657340457809111369f;   // = log(   [sqrt(2)]nearsmall  ) * F_LSS
	}

	// calc log
	sfXfrac=(sfX-1)/(sfX+1);
	sfHfsq=sfXfrac*sfXfrac;
	sfRet=K0*sfHfsq+K1;
	sfRet=sfRet*sfHfsq+K2;
	sfTemp=sfHfsq+S0;
	sfTemp=sfTemp*sfHfsq+S1;
	sfTemp=sfTemp*sfHfsq+S2;
	sfRet=sfRet/sfTemp;
	sfRet*=sfXfrac;


	if ( sfTemp2.st._LL != 0x0 ) {
		sfRet -= sfTemp2._F;
	}

	sfRet+=iK*F_LSS;




	return sfRet;
}

